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Abstract 



The model of laminated wave turbulence presented recently unites 
both types of turbulent wave systems - statistical wave turbulence (in- 
troduced by Kolmogorov and brought to the present form by numerous 
works of Zakharov and his scientific school since nineteen sixties) and 
discrete wave turbulence (developed in the works of Kartashova in 
nineteen nineties). The main new feature described by this model is 
the following: discrete effects do appear not only in the long-wave part 
of the spectral domain (corresponding to small wave numbers) but all 
through the spectra thus putting forth a novel problem - construction 
of fast algorithms for computations in integers of order 10 12 and more. 
In this paper we present a generic algorithm for polynomial disper- 
sion functions and illustrate it by application to gravity and planetary 
waves. 

Math. Classification: 76F02, 65H10, 65Yxx, 68Q25 
Key Words: Laminated wave turbulence, discrete wave systems, 
computations in integers, transcendental algebraic equations, com- 
plexity of algorithm 
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1 INTRODUCTION 



Statistical theory of wave turbulence begins with the pioneering paper [T| of 
Kolmogorov presenting the energy spectrum of turbulence as a function of 
vortex size and thus founding the field of mathematical analysis of turbu- 
lence. Kolmogorov regarded some inertial range of wave numbers between 
viscosity and dissipation, k < k < k\ for wave numbers k where k = \k\, and 
suggested that in this range turbulence is locally homogeneous and isotropic 
which, together with dimensional analysis, allowed Kolmogorov to deduce 
that energy distribution is proportional to k~ 5 ^ 3 . 

Kolmogorov's ideas were further applied by Zakharov for construction of 
wave kinetic equations [2] which are approximately equivalent to the initial 
nonlinear PDEs: 

A\ = J |K ( i23)| 2 5(Wi-CU 2 -CJ3)5(fc 1 -4-4)(^2v4 3 -AiA 2 -AiA 3 )dA : 2 dA : 3 

for 3-wave interactions, and similar equations for i-wave interactions where 8 
is the Dirac delta-function and Vm.i) is the vortex coefficient in the standard 
representation of nonlinearity in the initial PDE: 



s,W(fa + *» + - + *J jMl ..., 1< . (1) 

d\OJi + U0 2 + ... + UJi) 

The main idea of the wave turbulence theory is to take into account only 
resonant interactions of waves described by 



w(fci) ±u{k 2 ) ± ... ± u(k n+1 ) = 0, ^ 
k 1 ±k 2 ± ••• ± k n +i = 

Statistical wave turbulence theory deals with real solutions of Sys.(j2J) 

and one of its most important discoveries in the statistical wave turbulence 
theory are stationary exact solutions of the kinetic equations first found in |3] . 
These solutions have the form k~ a with a > and are now called Zakharov- 
Kolmogorov (ZK) energy spectra (see Fig.l). 

The left part of Fig.l, the so-called finite length effects, have been stud- 
ied in the papers of Kartashova [3] where properties of integer solutions of 
Sys. (PJ) have been studied. It was proven in particular that the spectral space 
of the discrete wave system is decomposed into the small disjoint groups of 
waves showing periodic energy fluctuations (nodes connected by lines at the 
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right panel of Fig. 2) and many waves with constant energy, depicted by 
blue diamonds. 

The most important result of the theory of laminated wave turbulence |Hj 
is following: discrete effects do appear not only in the long-wave part of the 
spectral space (corresponding to small wave numbers) but all along the wave 
spectra. From the computational point of view this theory gives rise to a 
completely novel problem: construction of fast algorithms for computations 
of integer solutions of Sys. (j2J) in integers of order 10 1 2 and more. For instance, 
for 4-wave interactions of 2-dimensional gravity waves this system has the 

form 

yki + \fk~2 = Vh + \fk\-, k 1 + k 2 = k ? , + fc 4 , 

where ki = (m^n,), Vi = 1,2,3,4 and k{ = \ki\ = \/mf + nf. It means that 
in a finite but big enough domain of wave numbers, say \m\, \n\ < D ~ 1000, 
direct approach leads to necessity to perform extensive (computational com- 
plexity D 8 ) computations with integers of the order of 10 12 . 

Importance of the discrete layer of laminated turbulence is emphasized 
by the fact that there exist many wave systems described only by discrete 
waves approach (for instance, wave systems with periodic or zero boundary 
conditions). A sketch of the first algorithm for the problems of laminated 
turbulence is given in 

In this paper we present generic algorithm for computing discrete layers 
of wave turbulent systems with dispersion function being a function of the 
modulus of the wave vector k, uj = uj{k). The main idea underlying our 
algorithm is the partition of the spectral space into disjoint classes of vectors 
which allows us to look for the solutions of Sys.(j2} in each class separately. 
In Sec. 2 we describe this construction in detail and with numerous examples 
because its brief description given in is often misunderstood by other 
researchers. In Sec. 3 the generic algorithm is presented with gravity waves 
taken as our main example while in Sec. 4 modification of this algorithm is 
given for the oceanic planetary waves. Results of the computations and brief 
discussion are given at the end. 
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2 DEFINITION of CLASSES 

For a given c £ N,c ^ 0,1,-1 consider the set of algebraic numbers R c = 
±/c 1//c , k E N. Any such number k c has a unique representation 

kc = iq 1/c ,l E Z 

where q is a product 

while pi,...p n are all different primes and the powers e±, ...e n E N are all 
smaller than c. 

Definition. The set of numbers from R c having the same q is called g-class 
Cl q (also called "class q"). The number q is called class index. For a number 
k{c) = 7<? 1//c , 7 is called the weight of fc( c ). 

The following two Lemmas are easily obtained from elementary properties 
of algebraic numbers. 

Lemma 1. For any two numbers ki,k 2 belonging to the same g-class, all 
their linear combinations 

ai&4 + a 2 k 2 , ai, a 2 E Z 
with integer coefficients belong to the same class q. 

Indeed, if k\ = 7ig 1 ^ c , k 2 = ^q 1 ^ then 

aiki + a 2 k 2 = ^iq\ /c + ^ 2 q\ /c = (71 + ^ 2 )q 1/c 

in other words, every class is a one-dimensional module over the ring of 
integers Z. 

Example 1. Let us take c = 2. Numbers k\ = and k 2 = y/lS belong 
to the same class 2: k\ = 2y / 2, k 2 = 3^/2. According to Lemma 1, their sum 
ki + k 2 belongs to the same class, and indeed, k\ + k 2 = \/50 = 5a/2. 
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Lemma 2. For any n numbers ki,k 2 ...k n belonging to pairwise different 
g-classes, the equation 

ki ± k 2 ... ± k n = 

has no nontrivial solutions. 

The statement of the Lemma 2 follows from some known properties of 
algebraic numbers and we are not going to present a detailed proof here. 
The general idea of the proof is very simple indeed: a linear combination of 
two different irrational numbers ^/gl, ^fq^ can not satisfy any equation with 
rational coefficients. For example, equation oa/3 + b\^E = has no solutions 
for arbitrary rational a and b. Indeed, taking the second power of both sides 
we immediately come to a contradiction, as the irrational number \/TE is not 
equal to any rational number. As for the general situation, one can apply 
the Besikovitch theorem jBj: 

Theorem. Let 

oi = hpi, a 2 = b 2 p 2 , a s = b s p s , 

where Pi,p 2 , p s are different primes, and bi, b 2 , b s are positive integers 
not divisible by any of these primes. If xi,x 2 , x s are positive real roots 
of the equations 

x ni - ai = 0, x™ 2 - a 2 = 0, x n ° - a s = 0, 

and P(x\, x 2 , x s ) is a polynomial with rational coefficients of degree less 
than or equal to n\ — 1 with respect to X\, less that or equal to n 2 — 1 with 
respect to x 2 , and so on, then P(xi, x 2 , x s ) can vanish only if all its 
coefficients vanish. 

Corollary. Any equation 

ai&i + a 2 k 2 ... + a n k n = 0, Oj G Z 

with ki, k 2 , ...k n belonging to pairwise different q-classes has no nontrivial 
solutions. 

1 /c 

Obviously, while h can be represented as 7^ , each number a{k{ = 
sgn(a i )\a i \>y i qy c . 
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Example 2. Let us take c = 2. Numbers k\ = and k 2 = VT2 belong 
to different classes: k\ = 2y/2, k 2 = 2\^3 and according to Lemma 2 the 
equation a\k\ + 02^2 = has only the trivial solution ai = a 2 = 0. 

The role of classes in the study of resonant wave interactions lies in the 
following theorem. 



Theorem 1. Consider the equation 

aiki + a 2 k 2 ... + a n k n = 0, Oj G Z 

where each ki = ^iq l J c belongs to some class a, G qi,q 2 ...qi, I < n. This 
equation is equivalent to a system 

91 ,l7 91 ,l + O gi ,27gi,2 + ••• + a gi,ni7gi,ni = 
^ a 92, 17^2,1 "I" a <?2,27<32,2 + ••• + a q2,n 2 lq2,ri2 = 

k a g ,,l7gi,l + a m,2l qi ,2 + ■■■ + a qi ,nil qi ,ni = 

Proof. Let us re-write the Equation grouping numbers belonging to same 
class as follows: 

( kqi, i J rk qi>2 +.. .-\-k qi >ni ) -\-(k q2j i-\-k q2>2 -\- . ..-\-kq 2jn2 )-\-...-\-{k qh i-\-k qh2 -\-...-\-k qhni ) = 

so that all numbers in the i-th bracket belong to the same class qi and all Oj 
are different. From Lemma 1 it follows that each sum k qu \ + k qi)2 + ... + k qi>ni is 
some number k qi of the same class qi and the equation can be re-written as 

k qi + k q2 + ---k^ = 

with all different qi. It immediately follows from Lemma 2 that k qi = 0, k q , 2 = 
+ ...k qi =0. 

Example 3. Now consider the equation 

01^ + 02^ + 03^ + 04^ + 05^ = (4) 
which is equivalent to 

201^ + 203^ + 303^ + 204^ + 405^ = (5) 
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and according to the Theorem 1 is equivalent to the system 

2eii + 3a 3 = 

2a 2 + 4a 5 = (6) 

CI4 = 

which is in every respect much simpler than the original equation. 

The computational aspect of this transformation is an especially impor- 
tant illustration to the main idea of the algorithm presented in this paper. 
Suppose we are to find all exact solutions of (jlj) in some finite domain 
1 < &i < D. Then the straightforward iteration algorithm needs 0(D 4 ) 
floating-point operations, = l..D,i = 1,2,3,4 (even ignoring difficulties 
with floating-point arithmetic precision for large D). On the other hand, 
solutions of Sys. (|5j) can, evidently, be found in 0(D) operations with integer 
numbers. 

3 EXAMPLE ONE: GRAVITY WAVES 

To show the power of the approach outlined above in practice, we proceed as 
follows. First we give a detailed description of the algorithm which is used 
to find all physically relevant four-wave interactions of the so-called grav- 
ity waves in a finite domain. We also estimate computational complexity 
and memory requirements for its implementation and present results of our 
computer simulations. In the next section we discuss reusability of this algo- 
rithm and transform it to solve a similar problem for three- wave interactions 
of planetary waves. Further we briefly discuss applicability of our algorithm 
to other wave-type interactions. 

3.1 Problem Setting 

The main object of our studies are four-tuples of gravity waves. A wave is 
characterized by its wave vector k, a two-dimensional vector with integer 
coordinates k = (m, n), m,n G Z, where m,n may not both be zero. We 
define the usual norm of this vector k = \k\ = V m 2 + n 2 and its dispersion 
function Wk = y/k. 
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Definition. Four wave vectors ki, k 2 , k%, k^ are called a resonantly interact- 
ing four-tuple if the following conditions are fulfilled: 

\fk[ + ^f¥ 2 = v^3 + Vh ^ 
k\ + k 2 = k 3 + k 4 



where k t = (mi,ni), Vi = 1,2,3,4 and ki = \k{\ = \/mf + nf. 

Sys. ((ZJ) is written in vector form and is equivalent to the following system 
of scalar equations: 

(m\ + nf) 1 ^ + (m\ + n^) 1 ^ = (m\ + n^) 1 ^ + (m\ + n^) 1 ^ 4 
m! + m 2 = m 3 + m 4 (8) 
ni + n 2 = n 3 + n 4 

Sometimes (especially in numerical examples) it is convenient to represent 
the solution four-tuple as 

Oil, n 1L )(m 2L , m 2L ) = (m 1R , n 1R )(m 2R , m 2R ) (9) 

We are going to find all resonantly interacting four-tuples with coordinates 
mi, rii such that —D < mi,rii < D, % — 1,2,3, 4 for some D G N. The set 
of numbers d G [— D, D] is further called the main domain or simply domain. 



3.2 Computational Preliminaries 
3.2.1 Strategy Choice 

Numerically solving irrational equations in whole numbers is always an in- 
tricate business. Basically, two approaches are widely used. 

The first approach is to get rid of irrationalities (for equations in radicals 
typically taking the expression to a higher power, re-grouping members etc.). 
For an equation like a^fx = b^fy this approach is reasonable: we simply raise 
both sides to power 2 and solve the equation a 2 x = b 2 y. (Some attention 
should be payed to the signs of a, b afterwards.) However, for Sys.(jEJ), con- 
taining four fourth-degree roots, this approach is out of question. 

The second approach is, to solve the equations using floating-point arith- 
metic, obtain (unavoidably) approximate solutions and develop some (do- 
main dependent) lower estimate for the deviation, which would enable us to 
sort out exact solutions with deviation due only to the floating-point. As an 
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example, consider the equation yfx = y in the domain < x,y < D. If x is 
not a square then \y/x — [y/x]\ > \/2^J\l/ D) — 1/8D, so each solution with 
smaller deviation is a perfect square. In other words, very small deviations 
are guaranteed to be an artefact of floating point arithmetic. 

This approach is more reasonable, though for Sys.(JBJ) the corresponding 
estimate would probably be not so easy to obtain. However, it has one cru- 
cial drawback, namely, its high computational complexity. Indeed, Sys.(JE} 
consists of 3 equations in 8 variables and exhaustive search takes at least 
0(D 5 ) operations, and many time consuming operations (like taking frac- 
tional powers) at that. 

Our primary goal is to find all solutions in the presently physically rel- 
evant domain D10 3 with a possibility of extension to larger domains. The 
algorithm should be generic, i.e. applicable to a wide class of wave types 
by simple transformations. Studying resonant interactions of other physi- 
cally important waves we may have to deal with even more variables, e.g. 
for inner waves in laminated fluid k = (m, n, I) and for four-wave resonant 
interactions the brute force algorithm described above has computational 
complexity 0(D 8 ). 

Clearly we need a crucially new algorithm to cope with the situation; and 
here classes come to our aid. 

3.2.2 Application of Classes 

Theorem 1, applied to the first equation of Sys.(JHJ) readily yields the following 
result. 

Theorem 2. Given an equation 

Vh = (10) 
with ti e N, U > 0, two situations are possible: 

Case 1: all the numbers ti, « = 1,2,3,4 belong to the same class Cl q . 
In this case Eq. (JlOJ) cab be rewritten as 

7iv^ + 72^g = 73^g + 74^g (11) 

with 71,72,73,74 £ N and q an i? 4 class index (i.e. a natural number not 
divisible by a fourth degree of any prime). 
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Case 2: the vectors belong to two different classes Cl qiJ Cl q2 . 

In this case Eq. (110)1 can be rewritten as 

71^ + 72^2 = 71^1+72^ (12) 
with 71, 72 G N and q±, qi being R4 class indexes. 

Now we concentrate on Case 1 being most interesting physically. For this 
case we can do computations class-by-class, i.e. for every relevant q we take 
all solutions of 71 + 72 = 73 + 74 such that jfq can be represented as a sum 
of squares jfq = m 2 + n 2 , \rii\ < D and for every decomposition into 

such sum of squares we check the linear condition (JB12). 

3.3 Algorithm Description 

At the beginning we have to compute a very important domain-dependent 
parameter we need for the computations. 

Notice that in the main domain — D < m,n < D every number under 
the radical t t < 2D 2 i.e, -ffq < 2D 2 . For a given q, i max (q) < {2D 2 /qf/ A . 

Definition. A number >-y max (q) is called class multiplicity and denoted Mul(q). 

For the main domain D = 10 3 , class multiplicities are reasonably small 
numbers, j max [l) = 37 being the largest. Class multiplicities for the majority 
of classes (starting with q = 125002) are equal to 1 - this fact will be later 
used to achieve a major shortcut in computation time. 

3.3.1 Step 1. Calculating Relevant Class Indexes 

Class indexes of the module R c as defined above are numbers not divisible 
by any prime in c-th degree, in our case c = 4 not divisible by 4-th power of 
any prime. We can further restrict relevant class indexes as follows. 

First, in (jlOj) every number under the radical ti = 'yfq must have a repre- 
sentation as a sum of two squares of integer numbers, ti = mf+n 2 . According 
to the well-known Euler's theorem an integer can be represented as a sum of 
two squares if and only if its prime factorization contains every prime factor 
p = Au + 3 in an even degree. As jf evidently contains every prime factor in 
an even degree, this condition must also hold for q. This can be formulated 
as follows: if q is divisible by a prime p = Au + 3, it should be divisible by 
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its square and should not be divisible by its cube. 



The implementation of this step is accomplished with a sieve-type pro- 
cedure. Create an array Ar q = [1, ...2D 2 ] of binary numbers, setting the all 
the elements of the array to 1. Make the first pass: for all primes p in the 
region 2 < p < \/2D 2 set to the elements of the array p 4 , 2p 4 , ...up 4 where 
k = [2D 2 /p 4 \. In the second pass, for all primes P4«+3 = 3 mod 4,p < 2D 2 
and integer factors a = l...a max such that ap < 2D 2 , do the following. If 
a ^ mod p then set the ap-th element of the array Ar q to 0. If a = 
mod p, then if a = mod p 2 then also set the ap-th element of the array 
Ar q to 0. 

Notice that in the second pass the first check should only be done for 
primes p < \f2D~ 2 and the second one - for p < V2D 2 . 

We create an array W q of "work indexes" . In the third pass, we fill it with 
indexes of the array Ar q for which the elements' values have not been set to 
in the first two passes. We also create an array of class multiplicities Mul q 
and fill it with corresponding class multiplicities (see previous subsection). 

Remark 1. All numbers q found above do have a representation as a sum 
of two squares; however, some do not have representation with \m\ < D and 
\n\ < D. However, we do not look for them now: they will be discarded 
automatically at further steps. 

The computational complexity of this step can be estimated in the follow- 
ing way. The number of primes < x is, asymptotically, tc(x) = x/\og(x), so 
their density around x is l/\og(x). The first pass takes [2D 2 /p 4 \ operations 
for each prime 2 < p < \J2D 2 so the overall number of operations can be 
estimated as 



As for the second pass, primes p = 4-u + 3 < 2D 2 constitute about a half of 
all primes and are evenly distributed among them. Sieving out by a prime p 
requires 0{2D 2 /p) + 0{2D 2 /p 2 ) + 0(2D 2 /p 3 ) = 0{2D 2 /p) operations which 
again boils down to overall 0(D 2 ) steps. 

Evidently, the third pass requires the same 0(D 2 ) operations. 

So the overall computational complexity of this step is 0(D 2 ). 
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Remark 2. It is not so easy to give a good estimate for the number of class 
indices ir d (D). Of course 0(D 2 /\n(D) < n cl (D) < 0(D 2 ) holds, and most 
probably ir c i(D) = 0(D 2 / \ia(D). (This is presently under study.) When- 
ever we need this number for estimating computational complexity of the 
algorithm, we presume vr c /(D) = 0(D 2 ) to be on the safe side of things. 
In our main computation domain D = 10 3 the number of class indices 
7r d (10 3 ) = 384145. 

3.3.2 Step 2. Finding Decompositions into Sum of Two Squares 

In 1908, G. Cornacchia ^Oj proposed an algorithm for solving the diophantine 
equation x 2 + dy 2 = 4p with p prime, p — Au + 1. This has been recently 
generalized to solving x 2 + dy 2 = m, m not necessarily prime [0]. To find 
all decompositions of a number 7 4 g into two squares we can use a simplified 
variant of this setting d — 1. A very efficient implementation of this algorithm 
can be obtained thanks to the following result 

Theorem 3. Let t 2 = — 1( mod m), < t < (m/2). Set r = m and r\ = t 
and construct the finite sequence {ri}, = qir i+ i + ri + 2, g« = |_?V r 'i+iJ , 
for < i < n — 1, where r > ri > ... > r n = 1 > r n+ \ = 0. If r\_ x > m> r\ 
then m = r 2 + r^ +1 . 

The proof is based on elementary number-theoretic considerations. Now 
it is evident that for each < t < (m/2) such that t 2 = — 1( mod m) we 
obtain one decomposition of m into two squares and the algorithm gives all 
decompositions with x > y. For our use, we also take symmetrical decom- 
positions x < y and also x — y if m — 2x 2 . The computational complexity 
T of the algorithm is, basically, the complexity of finding all square roots of 
— 1 modulo m and is logarithmic in m, i.e. T = 0(log(m)). 

Let Dec q be the maximal number of decompositions of 7 4 g, 7 = l...Mul(q) 
into sum of two squares. We create a three-dimensional array Ar£,[G q = 
l..Mul(q), D q = l..Dec q ,2] and for each G q store the list of decompositions 
( m G q ,D q ,nG q ,D q )- We also create a one- dimensional auxiliary array Ar DtoW 
storing the number of two-square decompositions for each weight. The num- 
ber of decompositions of an integer into sum of two squares can be estimated 
as 0(log(m)) using the following Euler theorem: 

Theorem. Let m be a positive integer, and let 

m = Tpl 1 ...pl k q[ 1 ...q t l l 
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be its factorization into prime numbers, where Pi = 1( mod 4) and = 3( 
mod 4). Then the number of essentially different decompositions of m into 
sum of 2 squares is equal to the integral part of 8/2 where 

i=i 

Now we see that filling the array Atd can be accomplished in 

T = 0(\og(q) + log(2 4 g) + ... + \og{Mul{qfq)) (13) 

steps. 

Using presentation ([TT!, Eq.(4.4.8.1)) 

n 

^ln(a£; + b) = n\na + \nT(b/a + n + 1) - lnT(6/a + 1) 

fc=i 

and the well-known 



//! ^ \/an(— ) n 



e 

we obtain T = 0(Mul(q)(logq + \og(Mul(q))) =0(Mul(q) log q). This is 
much less than 0(Mul(q) 3 ) (see the next subsection) and contribution of this 
step into the overall computational complexity of the algorithm is negligible. 



3.3.3 Step 3. Solving the Sum-of- Weights Equation 

Consider now the equation for the weights 

7i + 72 = 73 + 74 (14) 

with 1 < 7, < Mul(q) (see ITT|) . For convenience we change our notation to 
7il,72L (left) and 7i/j,72/? (right) and introduce weight sum S 7 = 71* + 72*. 
Without loss of generality we can suppose 

7lL < llR < 12R < 72L- (15) 

Notice that we may not assume strict inequalities because even for 7^ = 7^ 
there may exist two distinct vectors (mj,nj), (rrij,nj) with mf + nf = + 
n j = 7 4( ? either due to the possibility of representing 7 4 g as sum of two 
squares in multiple ways or even for a single two-square representation - to 
the possibility of taking different sign combinations (±|m|,±|n|) left and 
right. 



Now we may encounter the following four situations (see Fig.l a-d): 



15 



1- JlL < llR < 12R < 12L (Fig HI) 

The general, physically most interesting case. Every solution yields 
four waves with pairwise distinct modes. 

2. Til = Tifi < 12R = I2L (Fig© 

3. Til < Ti/? = 72L: < 72 l (Fig EI) 

4. Til = TiR = T2i? = T2L (FigCJ) 
The "most degenerate" case. 

The search is organized as follows. Each admissible sum of weights 
S 7 ; 2 < S 1 < Mul(q) is partitioned into sum of two numbers 5* 7 = 
Til + T2L, 1 < Til < T2L < Mul(q). Then the same number is par- 
titioned into sum of Ti-R;T2i?> Til < TiR < T2i? < T2L- Evidently, if 
Sj < Mul(q) + 1 then the minimal til is 1, otherwise it is 5 7 — Mul(q) 
(to provide Tii? < Mul(q)). The maximal til is always LS'- 7 /2j and similarly 
7ifl< LV2J- 

The computational complexity of this step can be estimated as T = 
0(Mul(q) 3 ) due to 0(Mul(q)) possibilities for each of three values S 7 , Til, Tifi- 
As Mul(q) = L(2D 2 /g) 1 /4j ; T = 0((2D 2 /g) 3 / 4 ). 

Remark. This step contains an evident redundancy. Indeed, the equation 
HUneed not be solved independently for each class. Instead, its solutions for 
all 5* 7 could be computed in advance and stored in a look-up table. However, 
this involves significant computational overhead (e.g. the lookup procedure 
includes computing the minimal Til = max(l, 5 7 — Mul(q)), which must be 
done for each class) wiping out the gains of this approach, at least for our 
basic domain D = 10 3 . Nevertheless, this approach should be kept in view 
if need for computations in much larger domains, say D = 10 6 , arises. 

Remark. The general case is not really so general - most classes have small 
multiplicities and then degenerate cases prevail. The overall distribution is 
given below: 



Case 


1 


2 


3 


4 


Classes 


24368 


57666 


13987 


63778 
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3.3.4 Step 4. Discarding "Lean" Classes 

In the main domain D < 10 3 we encounter 384145 classes. This sounds like 
a lot - however, most of these can be processed without computations or 
with very simple computations. Notice the simple fact that if a class has 
multiplicity 1, Sys.(JHJ) takes the form 

q = ml L + n \ L = m 2 L + n 2 L = m\ R + n\ R = m\ R + n\ R 
m 1L + m 2L = m 1R + m 2R (16) 
niL + n 2L = n lR + n 2R 

and for any nontrivial solution the four vectors (m^nj) should be pairwise 
distinct. In terms of the weight equation of the previous section it means 
that solutions, if any, have to belong to the fourth ("most degenerate") case. 
It is evident that no solution of Sys. fT6|) with pairwise distinct (m^n*) exist 
for q having few decompositions into sum of two squares: one (q = m 2 + m 2 ), 
two (q = m 2 + n 2 = n 2 + m 2 ) and three (q = m 2 + n 2 = n 2 + m 2 = I 2 + I 2 ). 
It can be shown by means of elementary algebra that this also holds for q 
having four decompositions. 

Remark. It is very probable that for classes of multiplicity 1 no nontriv- 
ial solutions exist, whatever the number of decompositions into sum of two 
squares. The question is presently under study. In the main domain D = 10 3 
we encounter 357183 classes of multiplicity 1 (1-classes). This is about 93% 
of all classes in the domain. Among them, the number of decompositions 
into sum of two squares is distributed as follows: 



Dec(q) 





1 


2 


3 


4 


5 


6 


7 


8 


Classes 


110562 


256 


138044 


163 


788 


86 


3 


8727 


2 


16595 


Dec(q) 


14 


16 


18 


20 


24 


26 


32 


9 


10 


12 




Classes 


38 


1015 


84 


1 


75 


1 


1 


31 


269 


2429 





Table 1. Distribution of decomposition number Mul(q) for 1-classes q in the 
main domain D = 1000 

It follows that 327911 1-classes can be discarded without any computa- 
tions at all and only 29272 must be checked for probable solutions. 

3.3.5 Step 5. Checking Linear Conditions: Symmetries and Signs 

Sum-of-weights equation solved and decompositions into sum of two squares 
found, we need only check the linear conditions to find all solutions. On the 
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face of it, the step is trivial, however some underwater obstacles have to be 
taken into account. Having found a four-tuple of vectors (ra», n$) satisfying 
the first equation of Sys.(JE} with both coordinates non-negative, solutions of 
the system will be found taking all combinations of signs satisfying 

±TOi ± m 2 = ±m 3 ± m 4 
±rii ± n 2 = ±n 3 ± n 4 

Even the straightforward approach does not need more than 2 8 comparisons, 
so this step does not consume very much computing time. However, a few 
points may not be overlooked in order to organize correct, exhaustive and 
efficient search: 

• For the system to represent a four-wave interaction, all the four waves 
must be pairwise unequal. For the first degenerate case of the weight 
diagram (FigJSJ) we must provide m±L ^ m± R and m 2R ^ m 2R . For the 
second one (FigEJ) - mm ^ m 2R and for total degeneration (FigEJ) - 
m 1L ^ m 2 L ^ m iR 7^ m 2R- 

• One and the same solution may not occur among the 256 sign combi- 
nations twice. First, this could happen due to some or being 
(evidently, the ± variation should not be done for any coordinate). 
Next, sign variation could lead to a transposition of wave vectors. For 
example, for q = 1 and 71 + 72 = 10 we obtain solutions 

{(0,-9)(0,49) (15,20)(-15,20) 
and 
(0,-9)(0,49) (-15,20)(15,20) 

which really represent one and the same four-tuple. 

• The set of solutions possesses some evident symmetries: if 

(m 1L ,n 1L )(m 2L ,m 2L ) =>> (m 1R ,ni R )(m 2R ,m 2R ) 
then, of course, 

{-m 1L ,n 1L )(-m 2L ,n 2L ) => (-m 1R ,n 1R )(-m 2R ,n 2R ) 

(m 1L ,-n 1L )(m 2L ,-n 2L ) => (m 1R , -n 1R )(m 2R , -n 2R ) 

(-m 1L ,-n 1L )(-m 2L ,-n 2L ) (-m 1R , -n 1R )(-m 2R , -n 2 R) 
Taking into account these points, an effective search is constructed easily. 
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4 EXAMPLE TWO: PLANETARY WAVES 



In this Section we demonstrate the flexibility of our algorithm. Namely, we 
solve essentially the same problem - finding all integer solutions in a finite 
domain —D < m,n < D for three-wave interactions and another wave type 
(the so-called "planetary waves" ). In this case c = — 2 and the main equation, 
corresponding to (jHJ), has the form 

1 

y/ m i+ n l (18) 

where \rrii\, |n$| < D. 

4.1 Steps that Stay 

• Step 1 - sieving out possible class bases - undergoes minimal changes. 
For c = —2, each q should be a square-free number and not divisible 
by any prime p = 4u + 3. Evidently, for this wave type the set of class 
indices is a subset of class indices of the previous section. 

• Step 2 - decomposition into two squares - can be preserved one-to-one. 
Indeed, there are sophisticated algorithms for representing square-free 
numbers as sums of two squares that are slightly more efficient than in 
the general case (one used in the previous section) but this step is not 
the bottleneck of the algorithm. 

4.2 Steps to be Modified 

• Step 3 - the weight equation is in this case 

1 1 1 

7i 72 73 

or 

7i72 

73 — 

7i + 72 

which has relatively few solutions in integers. Indeed, 
with multiplicity 1414 we obtain only 3945 solutions. 

Remark. For this example it makes sense to generate and store the 
set of triads (71, 72, 73) which constitute integer solutions of the Eq.fpHJj) 
for 1 < 7i < Mul(l) and for each class q just take its subset 1 < 7, < 
Mul(q). 



mi + ^2 = rriz 



(19) 
(20) 

even for class 1 
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• Step 4 - discarding "lean" classes - becomes trivial: no class with 
multiplicity 1 yields an integer solution of EH We need only consider 
63828 classes (g 638 28 = 499993, g 638 29 = 500009) from 243143 in the 
main domain. 

• Step 5 - checking linear conditions - is also much easier than in the 
previous example, i.e. one equation with three variables instead of two 
equations with four variables each. 

5 DISCUSSION 

Our algorithm has been implemented in VBA programming language; com- 
putation time (without disk output of solutions found) on a low-end PC (800 
MHz Pentium III, 512 MB RAM) is about 4.5 minutes for Example 1 and 
1.5 minutes for Example 2. Some overall numerical data for both examples 
is given in the Tables and Figures below: 



Domain 


< 200 


< 400 


< 600 


< 800 


< 1000 


Solutions 


263648 


800435 


932475 


1127375 


1389657 



Table 2. Example 1: Distribution of the number of solutions depending on 
the chosen main domain D. 

It is interesting that though the overall number of solutions grows sub- 
linearly as we extend the domain, the number of asymmetrical solutions 
(7i 7^ 72 7^ 73 7^ 74), physically most important ones, grows faster than 
linearly: 



Domain 


< 200 


< 400 


< 600 


< 800 


< 1000 


Solutions 


96 


344 


744 


1328 


2088 



Table 3. Example 1: Distribution of the number of the asymmetrical solu- 
tions depending on the chosen main domain D. 

Notice that considerable part of them (185 of the overall 2088) lie outside 
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of the D = 950 area, e.g.: 

(-150, -25) (990, 945) =>• (294, 49) (546, 871) 
where q = 37, 71 = 5, 72 = 15, 73 = 7, 74 = 13, 
(128,256)(990,180) (400, 200) (718, 236) 
where q = 20, 71 = 8, 72 = 15, 73 = 10, 74 = 13, 
(-80,-76)(980,931) (180, 171)(720, 684) 
where q = 761, 71 = 2, 72 = 7, 73 = 3, 74 = 6 

etc. As a whole, asymmetrical solutions are distributed not uniformly along 
the wave spectrum but are rather grouped around some specific wave num- 
bers. For instance, the first group of asymmetrical solutions (containing 8 
solutions) appears in the domain D = 50, with solution 

(_4, -4) (49, 49) (9,9)(36,36), 

and others, while in the domains D = 60, 70, 80, 90 there are no new asym- 
metrical solutions. The next new group (16 solutions) appears in the domain 
D = 100, and so on. From the physical point of view, asymmetrical solu- 
tions are the most interesting ones because they generate new wave lengths 
and, therefore, distribute energy through the scales. As it was pointed out 
quite recently [12], asymmetrical solutions play an extremely important role 
in wave turbulence. Indeed, no profound understanding of turbulence can be 
achieved without studying properties which is in our agenda. 



Numerical data for the case of planetary waves are given in the Table 4 
below: 



Domain 


< 200 


< 400 


< 600 


< 800 


< 1000 


Solutions 


1099 


3137 


5664 


8565 


11795 



Table 4. Example 2: Distribution of the number of solutions depending on 
the chosen main domain D. 



This data is presented graphically in Figs. IH1ITT1 Number of asymmet- 
ric solutions for Example 1 (gravity waves) and total solutions for Example 
2 (planetary waves) show smooth power growth and probably are asymp- 
totically power functions of the domain size D. On the contrary, the total 
solution number for Example 1 has an unexpected twist about D = 350 
(FigJHJ shown in detail at FigE]). This phenomenon is presently under study. 
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Notice that the algorithm presented here allows to find all solutions for 
wave vectors belonging to the same class. For three-wave interactions of 
arbitrary wave types this is always the case. For n-wave interactions with 
n > 3, however, interacting waves may belong to |_f J different classes [T3] . 
Consider for example the four-wave system 

Vh + Vh = Vh + Vh 
k 1 + k 2 = h + k A 

where k% and k^ belong to one class and k% and k^ - to another one, i.e. 
the first equation breaks up into two independent equations 

y/ki — \/^3 and \fk~2 = \fk&- 

In this case, a modified form of our generic algorithm can be applied. This 
will be dealt with in our next paper. 

Acknowledgement. E.K. acknowledges the support of the Austrian 
Science Foundation (FWF) under projects SFB F013/F1304. 
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Figure 1: Statistical Wave Turbulence Theory, \k\ > k$: Dependence 
of the energy E on the wave number k is presented in the right part of the 
picture; ZK-energy spectrum is shown symbolically. Left part of the picture 
is not explained by statistical wave turbulence theory. 
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n 



^ (m 2 + n 2 ) = k ^ k 



1/2 




m 



Figure 2: Discrete Wave Turbulence Theory, \k\ < k : 2D-spectral 
space is shown on the left panel, each node (n, m) of integer lattice cor- 
responds to the wave vector k = (n,m). Nodes corresponding to resonantly 
interacting waves are depicted by yellow squares and to non-interacting waves 
- by blue diamonds. On the right panel, dependence of the energy E on time 
t for both types of waves is shown 




Figure 3: Laminated Wave Turbulence Theory, arbitrary \k\: Discrete 
and statistical layers of turbulence co-exist in many wave systems. ZK- 
energy spectrum contains "holes" in the nodes of the integer lattice which 
are depicted by empty circles. 
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Figure 6: Weight diagram, Case 3: right degeneration 




Figure 7: Weight diagram, Case 4: total degeneration 
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Figure 8: Number of solutions in partial domains 




Figure 9: Twist range of Fig. |H] zoomed 
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Figure 10: Number of asymmetric solutions in partial domains 



1 20D0 




Figure 11: Example 2: Number of solutions in partial domains 
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